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Abstract. We consider the finite element solution of the vector Laplace equation on a 
domain in two dimensions. For various choices of boundary conditions, it is known that 
a mixed finite element method, in which the rotation of the solution is introduced as a 
second unknown, is advantageous, and appropriate choices of mixed finite element spaces 
lead to a stable, optimally convergent discretization. However, the theory that leads to 
these conclusions does not apply to the case of Dirichlet boundary conditions, in which both 
components of the solution vanish on the boundary. We show, by computational example, 
that indeed such mixed finite elements do not perform optimally in this case, and we analyze 
the suboptimal convergence that does occur. As we indicate, these results have implications 
for the solution of the biharmonic equation and of the Stokes equations using a mixed 
formulation involving the vorticity. 



1. Introduction 

We consider the vector Laplace equation (Hodge Laplace equation for 1-forms) on a two- 
dimensional domain Q. That is, given a vector field / on Q, we seek a vector field u such 
that 

(LI) curl rot — grad div It = / in Q. 

(Notations are detailed at the end of this introduction.) A weak formulation of a boundary 
value problem for this equation seeks the solution u in a subspace H C iJ(rot) fl if (div) 
satisfying 

(L2) (rot rot v) + (div It, div 1?) = (/, v), v E H. 

If H is taken to be if (rot) fl ii(div), the variational formulation implies the equation ( ILip 
together with the electric boundary conditions 

(L3) u-s = 0, divit = on dQ. 

Magnetic boundary conditions, u ■ n = 0, rottt = 0, result if instead the subspace H in the 
weak formulation is taken to be ii(rot) fl ii(div). (The terms electric and magnetic derive 
from the close relation of the Hodge Laplacian and Maxwell's equations.) If the domain 
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Q is simply-connected, both these boundary value problems are well-posed. (Otherwise, H 
contains a finite dimensional subspace consisting of vector fields which satisfy the boundary 
conditions and have vanishing rotation and divergence with dimension equal to the number 
of holes in the domain, and each problem can be rendered well-posed by replacing H with 
the orthogonal complement of this space.) 

Even when the domain is simply connected, finite element methods based on (11. 2p are prob- 
lematic. For example, on a non-convex polygon, approximations using continuous piecewise 
linear functions converge to a function different from the solution of the boundary value. 
See [21 § 2.3.2] for more details. A convergent finite element method can be obtained by 
discretizing a mixed formulation with a stable choice of elements. The mixed formulation 
for the electric boundary value problem seeks a G H^, u G i/(div) such that 

(cr, r) — (u, curl r) = 0, r G , 
(curl cr, v) + (div u, div v) = (/, v), v E H(div). 

On a simply connected domain, this problem has a unique solution for any vector field 
/; u solves (II. ip and (II. 3p and a = rotw. To discretize, we choose finite element spaces 
Tifi C H^, Vh C if(div), indexed by a sequence of positive numbers h tending to 0, and 
determine ct/i G S/^, it/i G Vh by 

(1.4) (o-fe,r) - ('a/j,curlr) = 0, r G E/j, 

(1.5) (curl cTfe, t;) + (div It,,, div v) = (/,«), v eVh- 

In order to obtain a stable numerical method, the finite element spaces E/^ and Vh must 
be chosen appropriately. A stable method is obtained by choosing E/j to be the Lagrange 
elements of any degree r > 1 and Vh to be the Raviart-Thomas elements of the same 
degree r (where the case r = 1 refers to the lowest order Raviart-Thomas elements). In 
the notation of [2], E^ x \4 = P^A" x V~A^ and the hypotheses required by |2] (the spaces 

belong to a subcomplex of the Hilbert complex if (div) with bounded cochain 

projections) are satisfied. From this it follows that the mixed finite element method is stable 
and convergent. Similar considerations apply to the magnetic boundary value problem, 
where the finite element spaces are E/, = E/j fl and Vh = VhH if (div) and the relevant 

Hilbert complex is if (div) L^. Another possible choice is to take E/, to be 

Lagrange elements of degree r > 1 and Vh to be Brezzi-Douglas-Marini elements of degree 
r - 1 (i.e., E/, X 14 = P^A° x Vr^iA^). This case is similar, and will not be discussed further 
here. 

We turn now to the main consideration of the current paper, which is the equation (11. ip 
with Dirichlet boundary conditions u = on dQ. This problem may of course be treated 
in the weak formulation (11. 2p with H = H^{Q] M^). In this case we may integrate by parts 
and rewrite the bilinear form in terms of the gradient (which, when applied to a vector, is 
matrix- valued) : 

(rot It, rot v) + (div u, div v) = (grad u, grad v), u,v E (f2; M^) . 
Thus the weak formulation (11. 2p is just 

(1.6) (gradu,gradt;) = (/,'y), v E H\n;R^), 
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for which the discretization using Lagrange or similar finite elements is completely standard. 

However, one might consider using a mixed method analogous to fll.4p -( lT3]) for the Dirich- 
let boundary value problem in the hope of getting a better approximation of cr = rot u, or 
when Dirichlet boundary conditions are imposed on part of the boundary and electric and/or 
magnetic boundary conditions are imposed on another part of the boundary. In fact, as we 
discuss in Sections S] and O a mixed approach to the vector Laplacian with Dirichlet bound- 
ary conditions is implicitly used in certain approaches to the solution of the Stokes equations 
which introduce the vorticity, and in certain mixed methods for the biharmonic equation. In 
the mixed formulation of the Dirichlet problem for the vector Laplacian, the vanishing of the 
normal component is an essential boundary condition, while the vanishing of the tangential 
component arises as a natural boundary condition. No boundary conditions are imposed on 
the variable cr. Thus, we define V/i = fl H{div), and seek ah & T,^, Uh ^ Vh satisfying 

(1.7) {ah, r) - {uh, curl t) = 0, t e T.h, 

(1.8) (cnilah^v) + {dwuh,dwv) = {f,v), v e Vh- 

Note that curl S/i ^ Vh, so there is no Hilbert complex available in this case, and the theory of 
[2] does not apply. This suggests that there may be difficulties with stability and convergence 
of the mixed method (ll.7l) - fll.8p . In the next section, we exhibit computational examples 
demonstrating that this pessimism is well founded. The convergence of the mixed method for 
the Dirichlet boundary value problem is severely suboptimal (while it is optimal for electric 
and magnetic boundary conditions). Thus, the difficulties arising from the loss of the Hilbert 
complex structure are real, not an artifact of the theory. 

However, the computations indicate that even for Dirichlet boundary conditions, the mixed 
method does converge, albeit in a suboptimal manner. While we do not recommend the 
mixed formulation for the Dirichlet problem, in Section E] we prove convergence at the sub- 
optimal rates that are observed and, in so doing, clarify the sources of the suboptimality. 
Theorem 13.11 summarizes the main results of our analysis, and the remainder of the section 
develops the tools needed to establish them. 

This analysis of the mixed finite element approximation of the vector Laplacian has impli- 
cations for the analysis of mixed methods for other important problems: for the biharmonic 
equation using the Ciarlet-Raviart mixed formulation, and for the Stokes equations using 
a mixed formulation involving the vorticity, velocity, and pressure, or, equivalently, using a 
stream function- vorticity formulation. As a simple consequence of our analysis of the vector 
Laplacian, we are able to analyze mixed methods for these problems, elucidating the sub- 
optimal rates of convergence observed for them, and establishing convergence at the rates 
that do occur. Some of the estimates we obtain are already known, while others improve 
on existing estimates. The biharmonic problem is addressed in Section H] and the Stokes 
equations in Section |3 

We end this introduction with a summary of the main notations used in the paper. For 
sufficiently smooth scalar-valued and vector-valued functions a and u, respectively, we use 
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the standard calculus operators 

, .da da , , ,da da, dui du2 du2 dui 
grada = ( — , — ), curl a = ( — , -— ), divn = — + — , rot w = — — . 

ox oy oy ox ox oy ox oy 

We use the standard Lebesgue and Sobolev spaces LPiVt), H\Q), Wp{Q), and also the spaces 
H{div,Q) and if (rot, ^2) consisting of vector fields u with divu in or rot u G L^, 
respectively. Since the domain Q will usually be clear from context, we will abbreviate these 
spaces as L^, H\ i7(div), etc. For vector- valued functions in a Lebesgue or Sobolev space, 
we may use notations like M^), although when there is little chance of confusion we will 

abbreviate this to simply HK The closure of C^{Vl) in if^, if(div), and if(rot), are denoted 
H^, i5f(div), if(rot). Note that if u G if(div), then the normal trace u ■ n E H~^^'^{dVL) 
and ii'(div) = { w G if(div) | u ■ n = on dVt }. Similarly, if(rot) = { it G if (rot) \ u ■ s = 
on dVt}. We write (-, ■) for the L^(f2) inner product (of either scalar- or vector-valued 
functions) and || ■ || for the corresponding norm. 

We shall also need the dual space of if(div), the space -ff(div)', normed by 

(1-9) ll^llH(div)' := sup J— . 

Clearly, 

(1.10) L^{VL] M^) C iJ(div)' C H-^{VL] M^) 

with continuous inclusions. 



2. Some numerical results 



We begin by considering the solution of the Hodge Laplacian (11. ip with electric boundary 
conditions (II. 3p using the mixed method (11.40 . (II. 5p . For the space S/i, we use Lagrange 
finite elements of degree r > 1 and for the space \4, Raviart-Thomas elements of degree r 
(consisting locally of certain polynomials of degree < r, including all those of degree < r — 1). 
These are stable elements and a complete analysis is given in [2]. Assuming that the solution 
is smooth, it follows from [21 Theorem 3.11] that the following rates of convergence, each 
optimal, hold: 

\\u-Uh\\ = 0{K-), II dw{u-Uh)\\ = Oih"-), \\a-ah\\ = 0{h'+'), \\ grad(a-aO|| = 0(/i^). 

Table 12.11 shows the results of a computation with r = 2. Note that the computed rates 
of convergence are precisely as expected. In the test problem displayed, the domain is 
Q = (0, 1) X (0, 1) and the exact solution is n = (cos yrxsin vr?/, 2 sinvrx cos vry). The meshes 
used for computation were obtained by dividing the square into 2" x 2" subsquares, n = 
1, 2, 4, . . . 128, and dividing each subsquare into two triangles with the positively sloped 
diagonal. Only the result for the four finest meshes are shown. Very similar results were 
obtained for the case of magnetic boundary conditions, and for a sequence of nonuniform 
meshes, and also for other values of r > 1. 

The situation in the case of Dirichlet boundary conditions is very different. In Table fI72\ 
we consider the problem with exact solution u = (sinvrxsinvri/, sinTrxsinvry). The finite 



MIXED FINITE ELEMENT APPROXIMATION OF THE VECTOR LAPLACIAN 



5 



Table 2.1. errors and convergence rates for degree 2 mixed finite element 
approximation of the vector Laplacian with electric boundary conditions. 



\\U - Uh\\ 


rate 


II div(tt-Mft)|| 


rate 




rate 


II curl(CT-a/0| 


rate 


2.14e-03 


1.99 


1.17e-02 


1.99 


2.16e-04 


3.03 


2.63e-02 


1.98 


5.37e-04 


1.99 


2.93e-03 


2.00 


2.70e-05 


3.00 


6.60e-03 


1.99 


1.34e-04 


2.00 


7.33e-04 


2.00 


3.37e-06 


3.00 


1.65e-03 


2.00 


3.36e-05 


2.00 


1.83e-04 


2.00 


4.16e-07 


3.02 


4.14e-04 


2.00 



element spaces are as for the computation of Table 12.11 except that the boundary condition 
of vanishing normal trace is imposed in the Raviart-Thomas space Vh- Note that the 
rate of convergence for a is not the optimal value of 3, but rather roughly 3/2. The rate 
of convergence of curia (i.e., the rate of convergence of a) is also suboptimal by roughly 
3/2: it converges only as h^^'^. For u, the convergence rate is the optimal 2, but the 
convergence rate for div u is suboptimal by 1/2. 

Table 2.2. errors and convergence rates for degree 2 mixed finite element 
approximation of the vector Laplacian with Dirichlet boundary conditions. 



\\U - Uh\\ 


rate 


II div{u-Uh)\\ 


rate 




rate 


II Curl(CT-cr/OII 


rate 


1.22e-03 


2.01 


1.55e-02 


1.58 


1.90e-02 


1.62 


2.53e+00 


0.63 


3.05e-04 


2.00 


5.33e-03 


1.54 


6.36e-03 


1.58 


1.68e+00 


0.60 


7.63e-05 


2.00 


1.85e-03 


1.52 


2.18e-03 


1.54 


1.14e+00 


0.56 


1.91e-05 


2.00 


6.49e-04 


1.51 


7.58e-04 


1.52 


7.89e-01 


0.53 



We have carried out similar computations for r = 3 and 4 and for nonuniform meshes and 
the results are all very similar: degradation of the rate of convergence by 3/2 for a and curl a, 
and by 1/2 for divn. However the case r = 1 is different. There we saw no degradation of 
convergence rates for uniform meshes, but for nonuniform meshes a converged in with 
rate suboptimal by 1 and curl a did not converge at all. 

The moral of this story is that the mixed finite element method using the standard elements 
is indeed strongly tied to the underlying Hilbert complex structure which is not present for 
the vector Laplacian with Dirichlet boundary conditions, and the method is not appropriate 
for this problem. Nonetheless the experiments suggest that the method does converge, albeit 
at a degraded rate. In the next section, we develop the theory needed to prove that this 
is indeed so, and also to indicate where the lack of Hilbert complex structure leads to the 
suboptimality of the method. 

3. Error analysis 

Theorem 13. H which is the primary result of this section, establishes convergence of the 
mixed method for the Dirichlet problem at the suboptimal rates observed in the previous 
section. In it we suppose that is a convex polygon endowed with a shape-regular and 
quasi-uniform family of triangulations of mesh size h. We continue to denote by E/^ C 
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and Vh C H{dw) the Lagrange and Raviart-Thomas finite element spaces of some fixed 
degree r > 1, respectively, with Vh = VhCi H{div). 

Theorem 3.1. Let u denote the solution of the vector Laplace equation (11. ip with Dirichlet 
boundary condition u = Q, and let a = rot u. There exist unique ah G ^h, Uh G Vh satisfying 
the mixed method (ll.7p -( |TT8|) . // the polynomial degree r > 2, then the following estimates 
hold for 2 < I < r (whenever the norms on the right hand side are finite): 

\\u — Uh\\ < Ch''\\u\\i, 

II div{u - Uh)\\ + ||a- - (T/JI + h\\ curl(cr - ah)\\ < C/i'"^/^(| ln/i|||it||,yi^ + ||u||i+i/2). 

If r = 1, the estimates are: 

W'u- — Uh\\ < C/i| In /ip(| In /i| 11^11 + 11^112), 

II dw{u — Uh)\\ + ||a" — cr/i|| + h\\ curl(cr — ah)\\ < Ch^^'^{\ In /i| ||tt||vi/^ + /i^/^||n||2). 

Note that the error estimate for u is optimal order (modulo the logarithm when r = 1), 
while (again modulo the logarithm), the estimate for divu is suboptimal by 1/2 order, and 
the estimates for a and curia are suboptimal by 3/2 order. This is as observed in the 
experiments reported above. Above and throughout, we use C to denote a generic constant 
independent of h, whose values may differ at different occurrences. 

The proof of this theorem is rather involved. Without the Hilbert complex structure, the 
numerical method is not only less accurate, but also harder to analyze. The analysis will 
proceed in several steps. First, in Section 13.21 we will establish the well-posedness of the 
continuous problem, not in the space x H[div), but rather using a larger space than 
with weaker norm for a. Next, in Section 13.31 we mimic the well-posedness proof on the 
discrete level to obtain stability of the discrete problem, but with a mesh- dependent norm 
on Hh- This norm is even weaker than the norm used for the continuous problem, which 
may be seen as the cause of the loss of accuracy. To continue the analysis, we then introduce 
projection operators into Vh and Hh and develop bounds and error estimates for them in 
Section 13.41 In Section 13.51 we combine these with the stability result to obtain basic error 
estimates for the scheme, and we improve the error estimate for Uh in Section 13.61 using 
duality. 

3.1. Preliminaries. First we recall two forms of the Poincare-Friedrichs inequality: 

(3.1) ||r|| < Cp||curlr||, t e H\ H^/^H < Cp|| grad^/^H, tfj e H\ 

Here denotes the subspace of functions in with zero mean. Similarly, we will use 
to denote the zero mean subspace of L^. 

Next we recall the Hodge decomposition. The space L^(f2; M^) admits a decomposition into 
the orthogonal closed subspaces cmlH^ and gradif^, or, alternatively, into the subspaces 
curl if ^ and gradif^. The decomposition of a given 1; e according to either of these 
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may be computed by solving appropriate boundary value problems. For example, we may 
compute the unique p G H"" and G such that 

(3.2) i> = curl p + grad 0, 

by a Dirichlet problem and a Neumann problem for the scalar Poisson equation, respectively: 

(curl p, curl r) = (v,curlr), r G n, 

(3.3) ^ ' ^ ' 

(grad 0, grad ^z^) = (i7, grad?/'), ip E H . 

Clearly, || grad0|| < \\v\\. If v G if(div), then satisfies the Neumann problem 



A(f) = div V in il, — = on d^l, / (f)dx = 0, 
on Jn 

so, by elliptic regularity, ||0||2 < C\\ divv|| if the domain is convex, and ||0||i < C\\ divv|| 
for any domain. 

We shall need analogous results on the discrete level. To this end, let Sh denote the space 
of piecewise polynomials of degree at most r — 1, with no imposed interelement continuity. 
Then the divergence operator maps onto Sh and also maps Vh onto S^, the codimension 
one subspace consisting of functions with mean value zero. The former pair of spaces is used 
to solve the Dirichlet problem for the Poisson equation, and the later is used to solve the 
Neumann problem. Each pair forms part of a short exact sequence: 

(3.4) 0^±h.^Vh^ Sh^O and ^ ±h ^ Vh ^ Sh ^ 0, 
respectively. 

The usual Raviart-Thomas approximate solution to Poisson equation A0 = g with Dirich- 
let boundary condition = is then: find Vh E Vh, (ph ^ Sh such that 

{vh,w) + {div w,(ph) = 0, weVh, {divvh,i') = {9,4^), ip e Sh- 

Define the operator grad^ : Sh — ?■ Vh by 

(grad^ 0, w) = -(0, divw), (I) E Sh, w E Vh.. 

From the stability of the mixed method, we obtain the discrete Poincare inequality ||0|| < 
Cpll grad^ 011, G Sh, with Cp independent of h. The solution {vh,(t>h) E Vh y< Sh of the 
mixed method may be characterized by 

(grad;, 0h, grad;, ^p) = -(g^^fj)^ ^jj e Sh, 

and Vh = gTadh(j)h- 

Corresponding to the first sequence in (13.40 . we have the discrete Hodge decomposition 

(3.5) Vh = curl Eh + grad^, Sh, 

and corresponding to the second, the alternate discrete Hodge decomposition 

(3.6) Vh = curl th + gradl Sh, 
where grad^ : Sh ^ Vh is defined by 

(grad^0,iy) = -{(f), divw), (p e Sh, we Vh- 
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Both of the discrete Hodge decompositions can be characterized by finite element compu- 
tations. For example, in analogy to fl3.3p . for given v G V/i we may compute the unique 
Ph G S/i and 0/i G Sh such that v = curl ph + grad^ (ph from the following finite element 
systems (one primal, one mixed): 

(curl p/j, curl r) = (f,curlr), r G S/i, 

(grad^ (f)h, grad^ ip) = {v, grad^ ip), ^ G Sh- 

3.2. Well-posedness of the continuous formulation. As a first step towards analyz- 
ing the mixed method, we obtain well-posedness of a mixed formulation of the continuous 
boundary value problem for the vector Laplacian. To do so, we need to introduce a larger 
space than for the scalar variable, namely 

S = {r G : curlr G i^(div)'}, 

with norm ||r||| = ||r|p + || curl r|||,^^.^^, (see (11. 9p ). The space S has appeared before in 
studies of the vorticity-velocity-pressure and stream function-vorticity formulations of the 
Stokes problem [10], and an equivalent space (at least for domains with C^'^ boundary) is 
used in The bilinear form for the mixed formulation is 

B{p, w; r, v) = (p, r) — (curl r, w) + (curl p, v) + (div w, div v), 

where (■, ■) denotes the pairing between if(div)' and H{dw) (or more generally between a 
Hilbert space and its dual.) Often, we will tacitly use the fact that if r is in H^, then 
(curl r, w) = (curl r, w). Clearly, 

\Bip,w;r,v)\ < 2(i|p||| + + 11^11^))'^'' P,rEE,w,vE i7(div), 

so B is bounded on (S x H{div)) x (S x if(div)). For r G S, we define Tq G by 

(curl To, curl?/') = (curl r, curl -i/^) , ip E n , 

Taking ip = tq shows that 

(3.7) II curl Toll < II curlr II ^(^i^y < ||r||E, r G S. 
It is also true that 

(3.8) ||r||s < C(||r|| + II curl roll), r G S. 
To see this, define G by 

(3.9) (0, div v) = (curlr, v) — (curl ro, v), v G if(div). 

This is well-defined, since divif(div) = L^, and, if diw vanishes, then v = curl?/; for some 
ip G H^, so the right-hand side vanishes as well. Clearly, 

(curlr,t;) = (curlro, v) + (</),divv) < (|| curlro|| + ||0||) ||f ||f/(div), v G i^(div). 

Choosing V G in ([33]) with divv = (p and ||v||i < C\\(p\\, we get \\(p\\ < C\\ curl(r-ro)||-i. 
This implies || curl rjlj^-^^j^y < C(||r|| + ||curlro||), thus establishing (13. 8p . We conclude 
from (13. 7p and (13.80 that the norm t t— ?■ ||r|| + || curlro|| is an equivalent norm on S. 
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Assuming that f E (or even H{dw)'), we now give a mixed variational formulation of 
the continuous problem. We seek a G S, u G i^(div), such that 

(cr, r) — (curl T,u) = 0, r G S, 
(curler,!?) + (div tt, div i;) = (/,v), v E H{div). 

We note that, if w G H^(il; M^) is the solution of the standard variational formulation (II. 6p 
and a = rot u, then a, u solve this mixed variational formulation. Indeed, u E H^{fl; M^) C 
H{div), a E L^, and, for v E H\n; M^) 

(curler, = (cr, rot?;) = (rot it, rot -y) = {f,v) — (div div t;). 

This implies that curler G H{div)', so cr G S, and, extending to G H{div) by density, that 
the second equation above holds. Finally 

(cr, r) = (rot u, r) = {u, curl r) 

for all r G L^, so the first equation holds. 

In the next theorem, we establish well-posedness of the mixed variational problem by 
proving the inf-sup condition for B, following the approach of [T]. Note that the theorem 
establishes well-posedness of the more general problem where the zero on the right hand side 
of the first equation is replaced by the linear functional {g,T), where g E S', and we allow 
/ G i7(div)'. 

Theorem 3.2. There exist constants c > 0, C < oo such that, for any (p, w) E T,x H(div), 
there exists (r, v) E x H[div) with 

(3.10) B{p, w; r, v) > c(||p||| + H||^(div)), 

(3.11) Iklls + ||f llH(div) < C'dlplls + \\w\\H{div))- 

Moreover, if w E curli^"^, then we may choose v E cmlH^. 

Proof. Define po G by (curl po, curl ?/)) = (curl p, curl ^/;), E H^. Next, use the Hodge 
decomposition to write w in the form w = curl/i + grad0, with /i G and (p ^ H^, and 
recall that 

(3.12) ||grad0|| < C|| divtw||. 
We then choose 

T = p — Sp, V = w + curl Po, 
where 5 is a constant to be chosen. Hence, 

B{p, W] r, v) = IIpII^ — 5{p, p) — (curl p, w) + 5(curl p, w) 
+ (curl p, w) + (curl p, curl po) + || div 
= IIpII^ + 6\\ curl /if - (5(p,p) + II curlpof + || dwwf. 

Recalling the constant Cp in the Poincare inequality (13. ip and choosing 6 sufficiently small, 
we obtain 

B{p,w;t,v) > ^||pf + (5 - (5^Cp/2)|| curl /if + || curlpof + || div wf 
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>C(||P||| + H||^(div)), 

where we have used the facts that = || curl/i|p + || grad^p, f l3.12p . and (13 .Sp in the 

last step. This estabhshes f lS.lOp . 

To estabhsh fl3.1ip . we observe that 

ll^'llH(div) < Hlk(div) + II curl Poll < ll^^lkcdiv) + IIpIIs 

by ([32]), while 

l|r|b < IIpIIe + sMj, < IIpIIs + SM, < c{\\ph + Nil), 

since ||p||i < C|| curl/i|| < C||it>||. 

To establish the final claim, we observe that if G curl if ^, then obviously v = w + 
curlpo £ curl if ^ □ 

Remark. Had we posed the weak formulation using the space x ii (div) instead of S x 
if(div), we would not have obtained a well-posed problem. 

3.3. Stability of the discrete formulation. In this section, we establish the stability of 
the mixed method (ll.7p - (ll.8p . guided by the arguments used for the continuous problem in 
the preceding subsection. Analogous to the norm on S, we begin by defining a norm on 
by Ikllv', = ||''"||^ + II curlrB , r G E/j, where 

II II (^''^) 
\\v\\^, := sup J — 7 . 

The bilinear form is bounded on the finite element spaces in this norm: 

\B{p,w;r,v)\ < 2(||p|||^ + ||^||^(di.))'/'(||r|||,^ + ||i;||^(div))'/', P,r e S,, w,v e Vh- 

For r G we define Tq G S/i by 

(curl To, curl ?/^) = (curl r, curl ?/^), ip G S/j. 

The discrete analogue of (13.70 again follows by choosing ip = curlro: 

II curlroll < II curlrll^^ < ||r||s^, r G S^. 

Next we establish discrete analogue of (13. 8p . that is, 

(3.13) ||r||s, <C(||r|| + II curlroll), r e E^. 

To see this, define G S*/! by 

(0, diw) = (curlr, v) — (curlro,!'), v G 

This is well-defined, since div Vh = Sh, and, if div v vanishes, then v = curl ip for some 
G S/i, so the right-hand side vanishes as well. It follows that || curlr||^/ < || curlro|| + ||0||. 

To bound ||0||, as in the continuous case, we choose v G with divv = (p and ||v||i < C||0||. 
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In the discrete case, we also introduce Hj^v, the canonical projection of v into the Raviart- 
Thomas space Vh (see f l3.22p ). so divil^i; = Pj^divu = and \\v — n^v\\ < Ch\\v\\i. 
Then 

||0f = (0,divi7,^^;) = (curlr,i7,^t;) - (curl tq, il,^^;) 

= (curlr, nj^v — v) + (curlr, v) — (curlro, n^v) 

< C/i(||r||i + II curlr||_i + || curlro||)||f ||i. 

Using the inverse inequality ||r||i < C/;.~^||r|| and the fact that ||i^||i < \\<P\\, gives the bound 
II0II < C'dlrll + II curl Toll), and implies (KUh . 

With this choice of norm, stability of the finite element approximation scheme is estab- 
lished by an argument precisely analogous to that used in the proof of Theorem 13.21 simply 
using the S/^ norm, the discrete gradient operator grad^, the discrete Hodge decomposition 
(13. 6p . the estimate (I3.13p . and the discrete Poincare inequahty, instead of their continuous 
counterparts. 

Theorem 3.3. There exists constants c > 0, C < oo, independent of h, such that, for any 
(p, it>) G S/j X Vh, there exists (r, v) G S/^ x Vh with 

B{p,w;r,v) > c(||p|||^ + \\w\\jj^^-^-^), 

||''"||Sh + ||''^||f/(div) <C(|| PIIsji ~l~ ||''^||_H'(div))- 

Moreover, if w ^ curlS/i , then we may choose v G curlS/j. 

Remark. Note that ||r||E,^ < ||r||s for r G T,h, but, in general equality does not hold. Had 
we used the S norm instead of the S/^ norm on the discrete level, we would not have been 
able to establish stability. 

3.4. Projectors. Our error analysis will be based on the approximation and orthogonality 
properties of certain projection operators into the finite element spaces: 

Ps, : ^ Sh^ Ps, : ^ P^^ : ^ S,, P^^ : H{dw) ^ Vh- 

For Ps^, we simply take the projection. By standard approximation theory, 

||s - Psh4lp < Ch^\\s\\wi, < / < r, 1 < p < oo. 

For and P±^, we use elliptic projections. Namely, for any r G H^, 

(curl Pe^t, curl p) = (curlr, curl p), p G S/,, (Ps^^r, 1) = (r, 1), 
and, for any r G 

(curl P|,^r, curl p) = (curlr, curl p), p G S/^. 
Then, by standard estimates, 

(3.14) ||c7-PE^a|| +/i||(7-Ps^a||i <C/i'||c7||i, 1 < / < r + 1. 
Moreover, 

(3.15) (curl[(T - Ps^a], v) < C/i||curl((T- Ps,^(t) II II div^ll, v e Vh, a e H\ 
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To prove this last estimate, we use the discrete Hodge decomposition f l3.5p to write v = 
curl7/i + grad^ iph, with •jh G and iph ^ Sh- As explained in § 13.11 the pair (grad^^ iph, i^h) G 
V/i X 5*^ is the mixed approximation of (grad?/', ■?/') where if) G solves A^/; = divi; in Vt. 
Since VL is convex, \\ip\\2 < C*|| divi;||. Therefore, 

(curl[a - PsftCr],'") = (curl[o- - Ps^o"], curl 7/, + grad^ = (curl[(T - PshC^], grad/^ ^/j) 
= (curl[(T - Ps^o-],grad;,?/'ft - grad^) < Ch\\ curl((T - Ps^o")!! ||?/;||2 
< C/i||curl(a-PE,a)||||divt;||. 

For P±^T, T G H^, we will use the estimate (due to Nitsche [15] for r > 2 and Rannacher 
and Scott [16] for r = 1; cf. also [3 Theorem 8.5.3]): 

(3.16) ||r - PtjWw^ < C/i'"ir||H.^, 1 < / < r + 1, 2 < p < 00, 
which holds with constant C independent of p as well as h. 

We define the fourth projection operator, P^^ : if (div) — )■ Vh, by the equations 
{Py^v, curlr + grad^ s) = (v, curlr) — (diw, s), r G S/i, s G Sh- 
in view of the discrete Hodge decomposition (13.61) . Py^^v G Vh is well defined for any v G 
i7(div). It may be characterized as well by the equations 

(3.17) (v -P^>^v, curlr) = 0, r G E/^, 

(3.18) (div['y - P^^v], s) = 0, sE Sh- 

Similar projectors have been used elsewhere, e.g., [TJ eq. (2.6)]. The properties of Py^ are 
summarized in the following theorem. 

Theorem 3.4. For v E H{div) and U E , 

(3.19) div P^^ V = Ps^ div v , P^^ curl U = curl P^^ U- 
Moreover, the following estimates hold 

(3.20) \\v - Py^v\\LP < Cph^\\v\\wi, 1 < / < r, 2 <p < cx), 

(3.21) \\div{v-P^^v)\\<Ch^\\dwv\\i, < / < r, 
whenever the norm on the right hand side is finite. 

Proof. The first commutativity property in (I3.19P is immediate from (13.181) . and the diver- 
gence estimate (I3.2ip follows immediately. For the second commutativity property, we note 
that curlP^^t/ G Vh and that, if we set v = curlt/ and replace Py^v by curlP^^?7, then the 
defining equations (I3.17p . (I3.18P are satisfied. 

To prove the estimate (I3.20p . we follow the proof of corresponding results for mixed 
finite element approximation of second order elliptic problems given in [11]. First, we in- 
troduce the canonical interpolant 77^ : if^(r2;]R^) — )■ Vh into the Raviart-Thomas space, 
defined through the degrees of freedom 

(3.22) V Jv-nwds, wEVr-i{e), '"''^'^^^ wEVr-2{T), 
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where e ranges over the edges of the mesh and T over the triangles. Then 

(3.23) \\v - nj^vllLP <Ch^\\v\\wi, 1 < / < r, 1 <p < cx), 
and, since div Uj^v = Ps^ divv, 

(3.24) \\dw{v - n)^v)\\LP < Ch^WdivvWwi, < / < r, 1 < p < oo. 

Writing v — Py^v = (v — Hj^v) + {Uj^v — P^^ v), it thus remains to bound the second term. 

From fl3.19p . div{Py^^v — TI^v) = 0, so Py^v — n)^v = curl ph for some ph G S^. Applying 

the decomposition (13. 2p . we have v — Uj^v = curl p + grad^/; for some p E and E H^. 
From f l3T7D . 

(curl p/i, curl r) = (curl p, curl r), r E S^. 
Thus, ph = P±^P and so satisfies the bound || curlp/iH^p < C|| curlp||ip given above in f l3.16p . 

Since 

(curl p, curl r) = {v — Uj^v, curl r) = (rot(t; — H^v), r), r E H^, 

p E satisfies — A p = rot('y — Hj^v). Using the elliptic regularity result of [IHl Corollary 
1], we have for 1 < p < oo that 



\p\\wi < Cp\\ Tot{v - v)\\w~^,p < Cp\\v -Uf^v 



LP- 



Following the proof of that result, the dependence of the constant Cp on p arises from the 
use of the Marcinkiewicz interpolation theorem for interpolating between a weak and an 
estimate. Using the explicit bound on the constant in this theorem found in pi Theorem 
VIII. 9. 2], it follows directly that Cp < Cp, where C is a constant independent of p. We 
remark that this regularity result requires the assumed convexity of Q, and does not hold for 
all 1 < p < oo if is only Lipschitz (c.f. [E]). Estimate fl3.20p follows by combining these 
results and applying (13. 23 p . □ 

Theorem 13.61 below gives one more property of Py^, inspired by an idea in [17]. To prove 
it we need a simple lemma. 

Lemma 3.5. Let p be a piecewise polynomial function with respect to some triangulation 
which is nonzero only on triangles meeting dVt. Then for any 1 < q < 2, 

\\p\\L.<Ch'/''~^"\\p\\L2, 

where the constant C depends only on the polynomial degree and the shape regularity of the 
triangulation. 

Proof. By scaling and equivalence of norms on a finite dimensional space, we have 

WpWlht) < C/i'/'^"ip|U2(T), p E Vr{T), 

where the constant C depends only on the polynomial degree r and the shape constant for 
the triangle T. Now, let 7^^ denote the set of triangles meeting dil. Then 

= IIPllL'3(r) - ^ IIpIIl2(t)- 
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Applying Holder's inequality we have 

and #7^^ < Ch~^ by the assumption of shape regularity. Combining these results gives the 
lemma. □ 

Theorem 3.6. Let 2 < p < oo. Then 

(3.25) {v - P^^v,cm\T) <Ch-'^/^'^/P\\v - Py^v\\Lp\\Tl r e ^h, veH{dw)r]LP. 

Proof. Define f G by taking the Lagrange degrees of freedom to be the same as those for 
r, except setting equal to zero those associated to vertices or edges in dQ. Then HfH < C||r|| 
and r — f is nonzero only on triangles meeting dQ. By (13.171) . 

{v — Py^v, curlr) = (v — Py^v, curl[r — f]). 

Let q = p/{p — 1), so 1 < q < 2. Applying Holder's inequality, the lemma, and an inverse 
inequality, we obtain 

{v - P,>^^t;,curl(r - f)) < \\v - P^^v\\lp\\ curl(r - t)\\li 

< C\\v - P^^v\\Lph^^^-^/P\\ curl(r - f)||i2 

<C\\v-P^^v\\Lr>h~'/^~'/P\\r-T\\L2, 

from which the result follows. □ 



3.5. Error estimates by an energy argument. Using the projection operators defined 
in the last subsection and the stability result of the preceding section, we now obtain a basic 
error estimate (which is not, however, of optimal order). 

Theorem 3.7. Let r > 1 denote the polynomial degree. There exists a constant C indepen- 
dent of the mesh size h and of p & [2, oo), for which 

Ik - (ThW + h\\a - (Jh\\i + \\u - Uh\\H{diw) 

< C W"^'^"^'" {pW'^Wn + M\i+in-i/v) > 2 < / < r, ifr>2, 
- |/,i/2-i/p (^p\\u\\^, + h^/^+^/i'\\u\\^ , ifr = 1, 

whenever the norms on the right hand side are finite. 

Proof. We divide the errors into the projection and the remainder: 

a -ah = {a - Pj^^a) + (Pj^^a - an), u-Uh = {u- Py^u) + {Py^u - Uh). 

Since, 

Ik - P^h(^\\ + ^Ik - ^s.c^lli < Ch'\\a\\t < Ch'\\u\\t+i, 1 < t < r + 1, 
and, by Theorem 13. 4[ 

\\u - Py^u\\H[diY) < Ch^\\u\\t+i, 1 <t <r, 
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the projection error satisfies the necessary bounds (without the p||it||vi/^ term on the right- 
hand side). 

Therefore, setting p = (Jh—PT.^'^ and w = Uh—P-^^u, it suffices to show that for 2 < p < oo, 
(3.26) IIpII + \\w\\H(di.)<C{\\a-P^^a\\+h\\a-P^^a\\i + h-^/'-^/P\\u-P^^u^^ 

Indeed, both cases of the theorem follow from (13.261) . Theorem 13. 4[ and the inverse inequality 
C^IIpIIi < IIpII- By the stability result of Theorem l3.3[ there exists (r, i>) G E/j x Vh satisfying 

B{p, W; T,V)>C {\\p\\l,^ + HllH(div)) . Iklls^ + ll*^llH(div) < C (IIpIIs,, + \\w\\H(div)) ■ 

By Galerkin orthogonality, 

B{p, w; r, v) = B{a - Pj^^a, u - Py^u; r, v) 

= (a- Ps^a, t) - {u- P^^u, curlr) + (curl(cr - Pj^^a), v), 

where we used the definition of B and (13.191) in the last step. Applying the Cauchy-Schwarz 
inequality. Theorem 13.61 and (I3.15p . we then obtain 

B{p, w- r, v) < C{\\a - Ps.af + h^\\ curl(a - P^,a)\\^ + h^'^-"^-"^^\\u - Pyu\\l,Y'' 

^{\\rr+\\v\\i,,„,y'\ 

Together, these imply (13.261) and so complete the proof of the theorem. □ 

Choosing p = | ln/i| in the theorem gives a limiting estimate. 
Corollary 3.8. The following estimates hold whenever the right hand side norm is finite: 

Ik - CFhW + h\\a - (T/Jli + \\U - Uh\\H{diw) 

'/i'-i/2(|ln/i|||w||^^^ + ||n||^+i/2), 2</<r, z/r > 2, 
h^'^ {\\nh\\\u\\w^^ + hy^\\uh) , ifr = l. 



< C 



For smooth solutions, choosing the maximum value of / = r in the corollary gives subop- 
timal approximation of a by order h^/"^^ and suboptimal approximation of u and divn by 
order h^^"^ (ignoring logarithms). In the next section, we show how to improve the error 
estimate for u to optimal order. The other estimates are essentially sharp, as demonstrated 
by the numerical experiments already presented. 

3.6. Improved estimates for u — Uh- Using duality, we can prove the following estimate 
for u — Uh in L^, which is of optimal order (modulo logarithms for r = 1). 

Theorem 3.9. These estimates hold whenever the right hand side norm is finite: 

_ ^ (h'\\u\\i, 2<l<r, ifr>2, 
""^ \h{\\nh\^l^uU.^ + \\uh), z/r = l. 
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Proof. Define G S, G if (div) by 

B{t, v; (j), w) = {v,u — Uh), r G S, v E H{dw). 

Thus w solves the Poisson equation — Aio = u — in Q with homogeneous Dirichlet 
boundary conditions, and (f) = — rot w. Under our assumption that f2 is a convex polygon, 
we know that w G H^, (f) G H^, and ||0||i + ||if II2 < C\\u — Uh\\. 

Choosing r = a — ah and v = u — Uh and then using Galerkin orthogonahty, we obtain 
\\u - Uhf = B{a -ah,u- uf, (f>, w) = B{a -ah,u- Uh] - PT.h<P, w - Py^w). 
The right hand side is the sum of following four terms: 

Ti = (o- - ah, (p - PT.,,(t>), T2 = -{u- Uh, curl[0 - Psfe^]), 

T3 = (curl[a - ah], w - Py^w), T4 = (div[w - Uh],(\.\v[w - Py^w]). 

We have replaced (■, ■) by the L^-inner products because G and a = rot it is in 
whenever the right hand side norm in the theorem is finite. For Ti, we use the Cauchy- 
Schwarz inequality, the bound ||0— -Psh0|| < C*^||0||i ^ C'^H'"— for the elliptic projection, 
and the estimate of Theorem 13.71 with p = 2 to obtain 

|Ti| < C J^'ll'""'"^ ^ ■"'^11' 2 < ^ < if r > 2, 
^~ yh{\\u\\i + h\\u\\2)\\u-Uh\\, ifr = l. 

Similar considerations give the same bound for T4. 

To bound T2, we split it as {Py^u — u, curl[0 — Pe^^]) and T2 = {uh — Py^u, curl[0 — Psh^])- 
The first term is clearly bounded by C/i'||w||i||'u — Uh\\, while, for the second, we use f l3.15p 
to find that 

\n\ < Ch\\ curl(0 - Ps, 0)1111 div(n, " Pvu)\\. 
Bounding div(n;j — Py^u) via Theorem 13.71 and (13.211) . we get 

\T2\ < Ch^\\u\\i\\u - Uh\\, 2<l<r, 

\T^\ < Ch{\\u\\i + h\\u\\2)\\u - Uh\\, r = 1. 

Finally, we bound T3. If r > 2, then we simply use the Cauchy-Schwarz inequality, the 
bound 

(3.27) \\w - Py^w\\ < Ch^\\w\\2 < Ch^\\u - Uh\\, 

and the p = 2 case of Theorem 13.71 to obtain 

IT^I <Ch^\\u\\i\\u-Uh\\, 2<l<r. 

If r = 1, then (13.271) does not hold. Instead we split T3 as (curl[cr — Ps,^cr],iu — P^^ w) + 
(curl[P2^0" — ah],w — P-^^w). Since \\w — P^^w\\ < Ch\\w\\i < Ch\\u — Uh\\, the first term 
is bounded by C/i||(t||i ||w — Uh\\ < (7/1111411211^ — Uh\\. For the second, we apply Theorem 13.61 
and (I3.20p to obtain 

|(curl[Ps, a - ah], w - P^w)\ < Ch-'/^~'/P\\w - P^w\\L4Pj,^a - ah\\ 

< Ch^/'^~^/Pp\\w\\w^]\Pj:,a - ah]\, 2 < p < 00. 
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By the Sobolev inequality, ||i(;||vi/^i < i^p||i(;||vi/2, where q = 2p/(2 + p) < 2. Moreover, from 

[T8] and a simple extension argument the constant Kp < Cp^^'^. Since ||iu||vv2 < C||iu||2 
with C depending only on the area of the domain, we obtain 

|(curl[Ps,a - a,], w - Py^w)\ < Ch^'^-^'Pp^'^P,:^a - - 2 < p < oo. 

By ( I3.14P and Theorem 13.71 with r = 1, 

||PE,a - ah\\ < \\a - Pj,^a\\ + \\a - a^W < C{h^'^~^'Pp\\u\\w^^ + h\\uh). 

Thus we obtain 

\T^\ < C (^h^~^/Pp^/^\\u\\w^ + /l3/2-VPp3/2||^||^j _ ^^11^ 2 < J9 < oo, 

and, by choosing p = \ \iih\ and noting that /i^/^l In/ip/^ is bounded, 

l^sl < Ch (I ln/i|^/^||it||vFi + ||n||2) \\u - Uh\\. 
The theorem follows easily from these estimates. □ 

4. The Ciarlet- Ravi art Mixed Method for the Biharmonic 

In this section, we show that the above analysis immediately gives estimates for the 
Ciarlet-Raviart mixed method for the biharmonic, including some new estimates which 
improve on those available in the literature. 

Given g G H~'^{Q) = {H'^{Q)y, the standard weak formulation of the Dirichlet problem 
for the biharmonic seeks U G such that 

{AU,AV) = {g,V), VeH\ 

Letting a := —AU G L^, we have Act = —g. Assuming that g G H~^{Q), as we henceforth 
shall, for Q a convex polygon, we have that U G H^{Q), a G H^{Q) and 

||f/||3+lkl|l<C||(7||_i. 

Hence (cr, U) e x satisfy 

(cr, r) - (curl U, curl r) = 0, t e H'^, 

(curia, curly) = (^,1/), V E H\ 

We note that a mixed formulation in these variables, but with spaces that are less regular, 
can also be given for this problem (c.f. [1]), but we shall not pursue this approach here. 

The Ciarlet-Raviart mixed method [6] for the approximation of the Dirichlet problem for 
the biharmonic equation using Lagrange elements of degree r, seeks G S^, Uh G T^h such 
that 

{ah, r) - (curl Uh, curl r) = 0, r G T.h, 
(curl(Tfe,curl\^) = (5(,\/), Veth- 
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This discretization has been analyzed in many papers under the assumption that f2 is a 
convex polygon. It is proved in [12] and [3] that for r > 2, 

\\U - Uhh < Ch'-WUWr+i, \W - < Ch'^-'WUWr+i. 

The former estimate is optimal, while the estimate for ||cr — cr/.J| is two orders suboptimal. 
The case r = 1 was analyzed in [17], where it was proven that 

\\U-Uh\\i < Ch^/^\lnh\^/^\\U\U, \W-ah\\ < Ch^^^\\nh\\\U\U. 

These estimates are suboptimal by 1/4 and 3/2 orders respectively (modulo logarithms) 
and require H'^ regularity of U. (As noted in [T7j, the same technique could be applied for 
r > 2 to obtain a 3/2 suboptimal estimate on \\a — (Th\\-) Below we improve the estimate 
on \\U — Uh\\i for r = 1 to an optimal order estimate (modulo logarithms), with decreased 
assumptions on the regularity of the solution U. 

We now show how to obtain all of these results from the analysis of the previous section, 
with only minor modifications. Let u = curl U. Then 

5((T,n;r, curly) = {g,V), {t,V) e x H\ 

Similarly, with Uh = cnilUh, 

B{ah, Uh- r, curly) = {g, V), (r, V) G S,, x t^. 

As above, set p = (yh — PT^^cr G w = Uh—Py^u G Vh- Note that w = curl fZ/j— curlP^^f/ G 
curl Eft. Subtracting the above equations and writing v for curly, we have 

B{p, w; r, v) = B{a - Pt.,,o-, u - Py^u; r, v), (r, v) G S?,, x curl S/,. 

Since the stability result of Theorem 13.31 holds over the space J^h x curlS/j, as stated in the 
last sentence of the theorem, we can argue exactly as in proof of Theorem 13.71 and conclude 
that the estimates proved in that theorem for the Hodge Laplacian hold as well in this 
context with one improvement. To estimate the term \\u — P^^u\\lp in fl3.26p . instead of 
using (13.201) . we note that \\u — P-^^u\\lp = \\ curl([/ — P±^U)\\lp and invoke (13.161) . In this 
way we avoid a factor of p. The improved estimates of Theorem 13.91 also translate to this 
problem, with essentially the same proof and a similar improvement. The dual problem is, 
of course, now taken to be: Find G S, lu G cmlH^ such that 

B(t, V] 0, w) = (v,u — Uh), T E T., V E curl H^. 

Thus w = curliy, where W solves the biharmonic problem A'^ W = rot(ii — Uh) E 
with Dirichlet boundary conditions, and (p = A W. The relevant regularity result, valid on 
a convex domain, is 

ll^lh + II0II1 < C\\Wh < C\\ rot(tt - Uh)\\-i < C\\u - Uh\\. 

The remainder of the proof goes through as before, with the simplification that now the 
terms T4 and T2 are zero, and the term \\w — P^^w\\lp can be bounded without introducing 
a factor of p as just described. The suppressed factors of p lead to fewer logarithms in the 
final result. Stating this result in terms of the original variable U instead of n = curlf/, we 
have the following theorem. 
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Theorem 4.1. Let U solve the Dirichlet problem for the biharmonic equation, a = —AU, 
and let Uh € S/i, Ch ^ denote the discrete solution obtained by the Ciarlet-Raviart mixed 
method with Lagrange elements of degree r > 1. If r > 2 and 2 < I < r, then the following 
estimates, requiring differing amounts of regularity, hold whenever the norms on the right 
hand side are finite: 



\a - ah\\ + h\\a - ah\\i < C 



h'-'\\U\W, 

h'^'^HmWwL-^ + \\U\\w/2) 
U-UH\\i<Ch'\\U\\i+i. 



Ifr = l, the estimates are: 

||c — o'hW + h\\a — (ThWi ^ C 



{\\uh + h\\uh), 

U-Uh\\i < Ch{\lnh\^/^\\U\\w^ + WUh). 



5. Stationary Stokes equations 

Another application in which the vector Laplacian with Dirichlet boundary conditions 
arises is the stationary Stokes equations, in which the vector field represents the velocity, 
subject to no-slip conditions on the boundary. A standard weak formulation (with viscosity 
equal to one) seeks u G H^{Q,M.^) and p E L"^ such that 

(gradtt, grad v) — {p,divv) = {f,v), v e H^{Q,,'R'^), 
(div u,q) = 0, g e L^. 

Mixed methods, such as we have discussed, have been used to approximate this problem, 
based on the vorticity-velocity-pressure formulation. For example, using the spaces defined 
in Section [3l the following weak formulation is discussed in [10]. Find cr G S, w G H{dw), 
p E Li^ such that 

(cr, r) — (curl r, w) = 0, r G E, 

(curler, v) — (p, diw) = {f,v), v G H{div). 

(divn, g) = 0, g G L^. 

This formulation is obtained just as for the vector Laplacian, by writing 

(grad u, grad v) = (rot u, rot v) + (div u, div v) 

and introducing the variable a = Totu. When / G L^(n;M^) and is a convex polygon, 
u G H'^{Q;M.'^), p G H^{Q), and a = Totu G H^{Q). Assuming this extra regularity, and 
setting u = curl U, and v = curl V, {a, U) G x satisfy the stream function- vorticity 
equations: 

{a, t) - (curl U, curl r) = 0, t e 
(curl a, curl V) = {f , curl V), V e H\ 
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Taking g = rot /, this formulation coincides with the mixed formulation of the biharmonic 
problem discussed in the previous section. 

We consider here the finite element approximation which seeks ah G S/i, Uh G Vh, Ph G Sh 
such that 

{ah, t) - {uh, curl r) = 0, r G S/j, 
(curlcr/j,v) - (ph, divv) = if,v), v G Vh- 
(divu, g) = 0, q e Sh- 

where the spaces T,h, Vh, and Sh are defined as above. The existence and uniqueness of the 
solution is easily established by standard methods. When / = 0, we get by choosing r = ah, 
V = Uh, q = Ph and adding the equations that ah = and div Uh = 0. Hence Uh = cmlUh, 
Uh G S/i, and choosing r = Uh, we see that cmlUh = 0. Since divVh = Sh, we also get 
Ph = 0. 

Error estimates for ||n — and ||cr — a/i|| are easily obtained by reducing the problem to 
its stream function-vorticity form and using the estimates obtained in the previous section. 
Letting Uh = cuilUh, and choosing v = curll^, V G S/i, we see that {ah, Uh) is the unique 
solution of the Ciarlet-Raviart formulation of the biharmonic with g = rot /. Hence, the 
estimates for a — a/i in Theorem 14. 1 1 remain unchanged, except that we can replace \\U\\s by 
In particular, we have the following theorem. 

Theorem 5.1. Let {u,p) solve the Dirichlet problem for the Stokes equation, a = rotn, and 
let Uh G Vh, ah G T^h, and ph G Sh denote the discrete solution obtained by the vorticity- 
velocity-pressure mixed method with r > 1 the polynomial degree. If r > 2 and 2 < I < r , 
then the following estimates, requiring differing amounts of regularity, hold whenever the 
norms on the right hand side are finite: 

\\a - ahW +h\\a-ahh<c[ ^ J"'j|';,||^^ + ||,||^^^/^) , 

\\u — Uh\\ < C/l'llwll;. 



If r = 1, the estimates are: 

— ah\\ + h\\a — ah\\i < C 



\\u\\i + /l||'u||2, 

h^/^{\\uU.^ + h^l^u\\,) 
u — Uh\\ < Ch{\ \Tih\^^'^\\u\\w^ + lli^lh)- 



The only item remaining is to derive error bounds for the approximation of the pressure. 
We obtain the following result, which gives error bounds that are suboptimal by 0(/i^/^). 



Theorem 5.2. If r >2 and 2 < I <r, then 
\\p-Ph\\<C 



h^-\\\u\\i + \\p\\i.^), 

h^'^''^ {WAwL + ll'"lb+l/2 + \\P\\l-l/2) 
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Ifr = 1, the estimates are 

\\p-Ph\\ < c 



\u\\i + h\\u\\2 H 

h}"{\\uU.^+h^l^uh+\\v\\l/2). 



Proof. From the variational formulation, we get the error equation 

{Ph - PshP^ div Vh) = {p- Ps^p, div Vh) + (curlfcr;, - a], Vh), Vh E Vh- 

We choose v G H^{^; M^) such that diw = ph — PsuP ^"^^ ^ C\\ph — Ps^p\\^ and take 
Vh = n'^v. We have that diw = div 77^ t; and ||i7^i;||//(div) < C'll'^'lli < C\\ph — PsuPWi so 

\\Ph - PsuPf = {Ph - PsnP, div nlv) = {p- Ps„p, div il^v) + {cm\[ah - a], nj^v), 
= {p- PshP,Ph - PshP) + (curl[(T/, - (y],nlv -v) + {ah - rot v) 
< C{\\p - PsM\ + h\\ cni\{ah - (J)\\ + \\ah - a\\)\\ph - PsMl- 
It easily follows that 

\\P-Ph\\ < C{\\p- Ps,M + \Wh - <y\\ + h\\ CUrl((Tft, - cr)||). 

The theorem follows directly by applying Theorem 15.11 and standard estimates for the error 
in the projection. □ 



A number of papers have been devoted to finite element approximation schemes of either 
the vorticity-velocity-pressure or stream-function-vorticity formulation of the Stokes prob- 
lem. In particular, the lowest order (r = 1) case of the method analyzed here was discussed 
in [9], (in which additional references can also be found). In the case of the magnetic bound- 
ary conditions, a = u ■ n = 0, the authors established stability and first-order convergence, 
which is optimal, for all variables. But for the no-slip boundary conditions u = 0, with 
which we are concerned and which arise much more commonly in Stokes flow, they observe 
in numerical experiments stability problems and reduced rates of convergence which are in 
agreement with the theory presented above. 

We close with a simple numerical example in the case r = 2 that demonstrates that 
the suboptimal convergence orders obtained above are sharp even for very smooth solutions. 
Our discretization of the vorticity-velocity-pressure mixed formulation of the Stokes problem 
then approximates the velocity u by the second lowest order Raviart-Thomas elements, the 
vorticity a by continuous piecewise quadratic functions, and the pressure p by discontinuous 
piecewise linear functions. We take Q to be the unit square and compute / corresponding to 
the polynomial solution velocity field u = (— 2x^(x — l)^y(2?/ — l){y — 1), 2y'^{y — l)^a;(2a; — 
l)(x — 1)), and pressure p = {x — 1/2)^ + iv — 1/2)^- The computations, summarized in 
Table \5A] indeed confirm the convergence rates established above, i.e., Uh converges with 
optimal order 2 to n in L^, while the approximations to a and curler are both suboptimal 
by 3/2 order and the approximation to the pressure p is suboptimal by 1/2 order. 
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Table 5.1. LP' errors and convergence rates for the mixed finite element ap- 
proximation of the Stokes problem for the vector Laplacian with boundary 
conditions ii-n = 0, tt-s = Oon the unit square. 



\\u-Uh\\ 


rate 


\\P-Ph\\ 


rate 




rate 


II curl(cr-(T,j)|| 


rate 


3.26e-04 


1.9 


2.34e-03 


1.3 


2.70e-03 


1.3 


1.67e-01 


0.2 


8.35e-05 


2.0 


8.05e-04 


1.5 


9.70e-04 


1.5 


1.24e-01 


0.4 


2.10e-05 


2.0 


2.74e-04 


1.6 


3.47e-04 


1.5 


8.96e-02 


0.5 


5.27e-06 


2.0 


9.39e-05 


1.6 


1.24e-04 


1.5 


6.42e-02 


0.5 
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